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Abstract. In this article we present our state of the art 
of fitting helioseismic p-mode spectra. We give a step by 
step recipe for fitting the spectra: statistics of the spectra 
both for spatially unresolved and resolved data, the use 
of Maximum Likelihood estimates, the statistics of the p- 
modc parameters, the use of Monte-Carlo simulation and 
the significance of fitted parameters. The recipe is applied 
to synthetic low-resolution data, similar to those of the 
LOI, using Monte-Carlo simulations. For such spatially 
resolved data, the statistics of the Fourier spectrum is as- 
sumed to be a multi-normal distribution; the statistics of 
the power spectrum is not a \ 2 with 2 degrees of freedom. 
Results for I = 1 shows that all parameters describing 
the p modes can be obtained without bias and with mini- 
mum variance provided that the leakage matrix is known. 
Systematic errors due to an imperfect knowledge of the 
leakage matrix arc derived for all the p-modc parameters. 
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1. Introduction 

In the past decade, helioseismology has been able to pro- 
vide the internal structure of the Sun and its dynamics. 
These inferences have been made possible by inverting 
the frequencies and rotational splitting of the pressure 
modes. The most commonly used technique for obtain- 
ing the p-mode parameters is to fit the p-mode spectra 
using Maximum Likelihood Estimators (MLE) assuming 
that the statistical distribution of the p modes in the 
power spectra is a \ 2 with 2 degrees of freedom (Woodard 
1984). The MLE with this statistics were first applied 
on helioseismology data by Duvall and Harvey (1986) 
and Anderson et al. (1990). This technique is used for 
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fitting spectra obtained with integrated sunlight instru- 
ments. For low- or high-resolution instruments, the (to, v) 
power spectra are commonly fitted assuming that each 
to spectrum has the same statistics as the for the inte- 
grated sunlight instruments (LOI instrument: (Appour- 
chaux et al. 1995; Rabello-Soares et al. 1997); GONG in- 
strument (Hill et al. 1996)). Unfortunately, none of these 
implementations are correct since the assumed statistics 
is wrong. Only Schou (1992) described a more correct way 
of fitting (to, v) diagrammes using not the power spectra 
but the complex Fourier spectra. 

The pioneering work of Schou (1992) has inspired this 
series of 3 articles for addressing our state of the art of 
fitting (to, v) diagrammes. In this paper (Part I), we de- 
scribe the statistics of the p modes, and how the MLE can 
be used in helioseismology. In Appourchaux et al. (1997) 
(hereafter Part II), we show how one can measure the 
mode leakage matrix and the noise correlation from the 
data which knowledge is required for using the Part I. In 
Appourchaux and Gizon (1998) (hereafter Part III), we 
applied these techniques to the LOI instrument of VIRGO 
on board SOHO (For a description of the performance of 
the instrument see Appourchaux et al. 1997). 

In this paper, we explain how the MLE can be used in 
helioseismology. In the first section, we recall the proper- 
ties of MLE. In the second section, we describe the statis- 
tics of the p-mode Fourier spectra. In this section, we have 
generalized the approach of Schou (1992), to any complex 
leakage matrices. We have also used complex matrices to 
generate the covariance matrices of the p modes and of 
the noise. In the third section we show how to use Monte- 
Carlo simulations for testing both the use of MLE and the 
model of the p-mode spectra, and then conclude. 

2. Maximum Likelihood Estimators 

Some of the properties of MLE were given by Toutain and 
Appourchaux (1994). We repeat them here for complete- 
ness. We also address 2 issues that were not covered in 
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their paper: are MLE biased?, and how significant are the 
estimated parameters. 

2.1. Fundamental properties 

The aim of this section is to introduce some definitions and 
properties of MLE. A comprehensive study of this area of 
statistics can be found, e.g. in Kendall and Stuart (1967). 
Given a random variable x with a probability distribution 
f(x,\), where A is a vector of p parameters. We define 
the logarithmic likelihood function I of TV independent 
measurements Xk of x as 



N 



]nL = t = -^2\nf(x k ,\). 



(1) 



fc=i 



where L is the likelihood. The main property of I is that 
the position of its minimum in the A-space gives an esti- 
mate of the most likely value of A, denoted hereafter as 
A. Hence A is the solution of the set of p simultaneous 
equations : 

dl 

— =0 with i = 1,2,..., p. (2) 



Moreover, in the limit of very large sample (N — > oo) 
this estimator A tends to have a multi-normal probability 
distribution. In this case, this estimator is asymptotically 
unbiased with minimum variance; which implies its expec- 
tation and variance are respectively: 



lim E(X) = A. 
lim er 2 (A) = ca 

N—>oo 



(3) 
(4) 



where Cu are the diagonal elements of the inverse of the 
Hessian matrix h, with elements: 



hij = E{ 



dt 



dXidX 



■)• 



(5) 



The covariances between any 2 components of A are given 
by the corresponding off-diagonal elements of the inverse 
matrix. Equation (5) is used when computing the so-called 
formal error bars on A; as a matter of fact according to the 
Cramer-Rao theorem, Eq. (5) gives only a lower bound 
to the error bars (Kendall and Stuart, 1967, reference 
therein). Toutain and Appourchaux (1994) showed that 
Eq. (5) is valid for most purpose in helioseismology. 

2.2. Biased or unbiased? 

The fact that MLE are asymptotically unbiased does not 
necessarily mean that this property is kept for a finite 
amount of data. As an example, it is well known that an 
estimator of the standard deviation (a) of N measurement 
of a normally distributed random variable x is given by: 



* 2 = 



1 



N 



N 



^2 



[xi — to) 



(6) 



where Xi is the i-th measurement of the random variable 
x and to is an estimate of the mean. It is well known that 
the a of Eq. (6) is unbiased. In this case, MLE would give 
the following estimator: 



N - 1 



® MLE — 



N 



(7) 



Clearly the MLE expression give a bias that vanish 
asymptotically for an infinite number of points. It is often 
difficult to derive explicit relation, similar to Eq. (7) be- 
tween the estimator and the finite number of data points. 
When analytical expression can not be found, we advice 
to use Monte-Carlo simulations to verify the unbiasness; 
an example for 1=1 splittings is given in Chang (1996) and 
Appourchaux et al. (1997). 

In any case MLE are intrinsically biased estimators be- 
cause they are also minimum variance estimators (Kendall 
and Stuart, 1967). It may be useful to find other estima- 
tors that do not bias the estimates (Quenouille, 1956); 
they might not necessarily have minimum variance. These 
estimators are yet to be found. 

2.3. Significance of fitted parameters 

When one uses Least Square for fitting data, one can test 
the significance of its fitted parameters using the so-called 
R test (Frieden, 1983). For MLE, a useful test can be used: 
the likelihood ratio test. It was first used by Appourchaux 
et al. (1994). This method requires to maximize the likeli- 
hood e~ e ^ of a given event where p parameters are used 
to described the line profile. If one wants to describe the 
same event with n additional parameters, the likelihood 
e -i{si p+n ) w jyj j^yg ^0 be maximized. The likelihood ra- 
tio test consists in making the ratio of the two likelihood 
(Brownlee, 1965). Using the logarithmic likelihood, we can 
define the ratio A as: 



ln(A) = l(n p+n ) - £{w p ) 



(8) 



If A is close to 1, it means that there is no improvement 
in the maximized likelihood and that the additional pa- 
rameters are not significant. On the other hand, if A <C 1, 
it means that £(Q p+n ) <C £(u) p ) and that the additional 
parameters are very significant. Wilks (1938) showed that 
for large sample size the distribution of — 21nA tends to 
the X 2 (n) distribution. 

3. The statistics of p-mode spectra 

3.1. Single mode 

It is well known that p modes are stochastically excited 
oscillators (Kumar et al. , 1988). The source of excitation 
lies in the many granules covering the Sun. The modes are 
assumed to be independently excited provided that their 
spatial scale is larger than the granule size (Chang, 1996). 
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From the equation of an oscillator, the statistics of the 
p-mode profile can be derived as: 



Cp 1 x dx 

^ + 2n 7 - + (2n)^ x = F { t) 



(9) 



where t is the time, x is the displacement, 7 is the damping 
term or the linewidth, v a is the frequency of the mode and 
F(t) is the forcing function. From this equation the Fourier 
transform of x can be written as: 



x(u) 



F{y) 



i'yu) 



(10) 



where x{y) and F(v) are the Fourier transform of x(t) 
and F(t). From the large number of granules, it can be 
derived that the forcing function is normally distributed. 
Therefore the 2 components (the real and imaginary parts) 
of the Fourier transform of the forcing function are also 
normally distributed. For the p modes, each component of 
the Fourier transform is normally distributed with a mean 
of zero, and a variance given by: 



1 



2(27t) 4 [K 2 -z/2) 2 + z/ 2 7 2 



(11) 



The square of the modulus of x{v), or power spectrum, 
has a x 2 with 2 degree of freedom statistics and its mean 
is given by Eq. (11). This is the p-mode profile that is 
usually approximated by a Lorentzian profile. Similarly 
other effects such as asymmetry can be introduced in the 
profile of Eq. (11). 

3.2. Unresolved observations 

Instruments integrating over the solar surface the velocity 
or the intensity signal observe a superposition of various 
modes of different degrees. They are mainly sensitive to 
the low-degree modes (I < 4). For a given I, they detect 
a mixing of azimuthal order m for which a visibility is 
prescribed (Toutain and Gouttebroze, 1994; Christensen- 
Dalsgaard and Gough, 1982). Most often they can only 
detect modes for which I + m is even. Since the Fourier 
components of the observed time series have a normal dis- 
tribution, and since the different m are uncorrelated, the 
statistics of the power spectra of unresolved observation 
is a x 2 with 2 degrees of freedom. Toutain and Appour- 
chaux (1994) gave an analysis of the problem associated 
with these observations; we will not repeat it here. 

3.3. Resolved observations 

When the solar image is resolved many more degrees can 
be detected making the data analysis somewhat more 
complicated. In order to extract a single l,m mode from 
resolved observations, one has to apply a specific spatial 
filter or weight to the velocity or intensity images. Most 
often these weights are such that imperfect isolation of the 



I, m mode is achieved; especially because the most com- 
monly used filters (spherical harmonics) are not orthogo- 
nal over half a sphere. This leads to the existence of other 
modes in the Fourier spectrum generated for a given /, m 
filter. Therefore, the observed Fourier spectrum is a lin- 
ear combination of the modes to be detected. This linear 
combination of the modes can be understood as modes 
leaking into each other spectrum: this is represented by 
the so-called leakage matrix. These leakages will produce 
correlations between the different Fourier spectra. These 
correlations will modify the statistics of the Fourier spec- 
tra, such that their power spectra cannot be described as 
a x 2 with 2 degrees of freedom. Therefore the statistics of 
the 21 + 1 power spectra of a given / cannot be derived from 
the product of 21+1 x 2 with 2 degrees of freedom as in Ap- 
pourchaux et al. (1995). Nevertheless, the real and imag- 
inary parts of the Fourier spectra will still be normally 
distributed; in other words, the Fourier spectra have a 
multi-normal distribution defined by a covariance matrix. 
This fact will be used to derive the statistics of the obser- 
vation. The covariance matrix is the sum of the noise and 
mode covariance matrices, which are not necessarily the 
same. Last but not least, the theoretical probability distri- 
bution has to be generated using the previous covariance 
matrix. 

In summary, to understand the statistics of resolved 
observation, one has to follow four steps: 

— Compute leakage matrices, 

— Compute mode covariance matrices (related to the 
leakage) , 

— Compute noise covariance matrices, 

— Generate the likelihood from the theoretical probabil- 
ity distribution 

Each step is described in detail hereafter. 



3.3.1. Leakage matrices 

Due to the spherical symmetry of the Sun, the most likely 
weights to be used to isolate the modes are the spheri- 
cal harmonics Y^ m . Here we generalize the approach to 
any weight W; >m . The result is the observation of Fourier 
spectra y l m (y) that are related to what we want to de- 
tect, i.e. the Fourier spectra of the individual Fourier spec- 
trum x l m ,{v), by the so-called leakage matrix (Schou 1992, 
Schou and Brown 1994). The following expression can be 
derived for as many different degrees as needed; for sim- 
plicity we wrote it for 2 different degrees I, I' as: 



= C(U') 



x 



(12) 



where x(u) and y{v) are 2 complex vectors made each of 
21 + 21' + 2 component: 21 + 1 components for I, 21' + 1 
components for I' and C^ 1 ' 1 ' is the leakage matrix of both / 
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and /'. The dimension of C W is (21 + 21' + 2) x (21+21' + 2). 
The coefficient of the leakage matrix can be computed as: 



7(1,1') 



,0,0 

m.m' 



(13) 



with: 



ix) = [ . 

n.m' / 
JV 



Y Vim , (6, <j>)A(e,<f>) sin 6d0d<j> (14) 



where m = —I, ■ ■ .,1, m' = —I', . . . , I', * denotes the com- 
plex conjugate, 6, <p are the angles in a spherical coordi- 
nate system, T> is the integration domain, jTO / is the 
generalized velocity or intensity perturbation of the mode 
(l',m'), A(6,4>) is an apodization function, n; iTO is a sen- 
sitivity correction factor associated with W\_ m . The ratio 

ensures that Cm,m = 1 • The apodization function A is the 
product of 3 different function as: 



A(6,cb) = A n (6,cf>)A d (0,cj>)A a (6,c!>) 



(15) 



for the GONG instrument (Hill, 1997, private communi- 
cation). In both cases, this is not produced by the obser- 
vation techniques but by the data analysis techniques. 

If the weight functions W^ m and the observed pertur- 
bations Y; jTO have the same symmetry properties as the 
spherical harmonics Yi tTn (or if Wi >m — Yi_ m — Y;, TO ), the 
leakage matrix is real as shown by Schou (1992). In addi- 
tion the leakage elements of C^f^, are zero if l+m+l'+m' 
is odd; this is the case when the Sun is not tilted with re- 
spect to the observer's axis of reference (P=0, -8=0). If 
the axes of reference of Wi. m differ from that of the Y/ iTO 
these 2 properties can be lost. For instance, an incorrect 
orientation of the Sun axis with respect to the detector 
axis could lead to a complex leakage matrix; or a Sun 
seen at an angle B ^ give a real leakage matrix with 
non-zero elements with l + m + l' + m' odd. This latter 
property has been used by Gizon et al. 1997 to infer the 
inclination of the Sun's core. 

Equation (13) is valid when the size of the pixel is small 
compared with the spatial scale of the degree. When the 
pixels are larger, one should write the following: 



A n is the natural apodization function due to the way 
the images are obtained: for intensity this is the limb 
darkening (1(h)), and for velocity the projection factor 
(fi = sin 9 cos (f>). Ad is the data analysis apodization: for 
data re-mapped on the Sun'surface it is unity; for no re- 
mapping, it is the projection factor (^ = sin # cos <?!>). A a 
is the artificial apodization that can take into account the 
non-linear velocity (or intensity) response of the instru- 
ment over the solar disk, or can help to reduce limb ef- 
fects. Here we must point out that the leakage matrix has 
a useful property such as: 



>0,O 



,0',0* b 



0,0* d 1 ' 1 ') 

m,m u m,m' 



m'm' m'm 



(16) 



It shows that C^ 1 ' 1 ' is in general not hermitian nor sym- 
metrical. Nevertheless, when Wi. m — Y^ m , it is possible 
with a proper sensitivity factor correction of W^ m to have 
such a property. In this case the sensitivity correction is 
given by: 



ni, m = ]jf v Y£ m <f>)Yi,m (*> <t>) A i e , <t>) ™ 9dBd4 (17) 

which is the 'natural' normalization factor of the pertur- 
bation Y „, . Of course in this latter case, we have: 



C (U') =c o'.o* 

m.m' m'm 



(18) 



Unfortunately, the leakage matrix does not always have 
such a nice property, especially because Wi. m ^ Yi ;m . 
This was the case for the ground-based Luminosity Os- 
cillations Imager (LOI) (Appourchaux et al. , 1994) and 



E(Z,m)* , m') 
i W l Vi 

°m,m' „ ^ U',m')*~(l',m') 

where the yn are given by: 



(19) 



yf' m,) =l YF'(Q, <j>)A(6,<P) sin 6ded<j> (20) 



where T>i is the area defined by the i-th pixel and wf '"^ 
is the weight applied to the i-th. pixel to extract the /, m 
mode. Equation (19) is the more general form used for the 
LOI (Appourchaux and Andersen, 1990). As a starting 

point, the wf ,m ^ can also be taken as the y^ ,m \ 
3.3.2. p-mode covariance matrix 

To compute the covariance of the complex vector y(v) as 
a real number we form the vector z y (v) defined as: 

4(z,) = (Rc(y T ),Im(y T )) 

In absence of noise, the covariance matrix M(i/) of the 
vector z y (v) can be generated using a complex notation: 



K ' \-Mi(v) M r (v) J 



(21) 



M^'' ) is a super matrix where M r (v) and Mi(v) are 
the real and imaginary parts of a complex matrix M^ 1 ' 1 ) 
which elements are given by: 

Mt i >)= e e ee^w (22) 

I"— 1,1' m" — — l" 
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where f l m n(v) is thc variance of the l",m" mode which 
profile is given by Eq. (11), in which vq is a function of m. 
The real and imaginary parts of Eq. (22) will give respec- 
tively the covariance of the real (or imaginary) part of y, 
and the covariance between the real and imaginary part 
of y. It is obvious from Eq. (22) that MS 1,1 is hermitian. 

Schou (1992) gave an equation similar to Eq. (22) for 
a real leakage matrix and for a single degree. Here we add 
a subtlety to the formulation of Schou (1992), the matrix 
MW)(i/) can be decomposed as follows: 



v ' \—\n{u) v(z/) J \yj 



-w T (z/) 



» 



(23) 



where T is the transpose of a matrix. The elements of v 
and w are given by: 



ft(^^>)<l) 



(24) 



(25) 



We will see later on that this decomposition is of prime 
importance for understanding the statistics of the obser- 
vation. 

3.3.3. Noise covariance matrix 

Unfortunately, the observed vector y(v) include a noise 
contribution. Due to the way the data are combined, the 
noises between thc different 21 + 21' + 2 components of this 
vector are also correlated. Schou (1992) gave the correla- 
tion matrix when thc filter used are spherical harmonics 
Yi m . A more general formulation can be written as: 



-Bi{u) B r {v) 



to know a model of the solar noise. An easier way to under- 
stand the noise correlation is to built the ratio covariance 
matrix or 'pseudo' noise leakage matrix 1Z as: 



(i,0 



B 



(i,0 



B (i',0 

m' ,m' 



(29) 



Here we can see the similarity between TZ and C. In veloc- 
ity, the granulation noise is rather low at the center of thc 
disk and then increases towards the limb; the meso- and 
super-granulation exhibits somewhat different or comple- 
mentary center-to-limb variations. In intensity, thc gran- 
ulation noise is a function of thc number of granules; thc 
noise is larger at the center of the disk and decreases slowly 
towards the limb. In addition the solar noise in inten- 
sity has no contribution from mesogranulation (Frohlich 
et al. , 1997), making the spatial dependence of the noise 
almost independent of frequency across thc p-mode range. 
This is not the case in velocity where mesogranulation still 
contributes to the noise in the p-mode range. Therefore in 
intensity the apodization a is closer to A than in velocity, 
making the ratio covariance matrix TZ^ 1 ' 1 ) very close to the 
leakage matrix C^ 1,1 \ Although 7Z^ 1,1 ** is not mathemati- 
cally useful, it is a matrix easy to visualize and understand 
(See Part II) . The ratio matrix has some properties of the 
leakage matrix like being not necessarily hermitian. This 
is not the case of B^ 1,1 ' which is hermitian by definition. 

Again, when the size of the pixel is large compared 
with the spatial scale of the degree, Eq. (27) is rewritten 
as follows: 



b (u '\m = Y 



w. 



■■(1,1 



2) (l',m') u I \ 

'w\ 'bi(y) 



(30) 



where hi is the variance of the noise of pixel i. Equation 
(26) (30) is the more general form used for the LOI. 



B^' J ) is a super matrix where B r {v) and Bi(y) are the 
real and imaginary parts of the complex matrix B^ 1,1 \ 
The dimension of fiW) is (2/ + 2V + 2) x (21 + 21' + 2). Its 
elements are given by: 



m 

with 



l '\ = [ 

.m' I 

JV 



wt tm (0,4>)w v , m ,{6,<p) 



nv, r . 



a(6,<j>) sin 0d0# (27) 



(28) 



where a is an apodization function which characterizes 
through CTq(0, <j), v) how the noise varies over the solar im- 
age, assuming that the noise is uncorrelated between dif- 
ferent points on thc Sun; A a ,A^ are defined in Eq. (15). 
When the instrumental noise is low, a is derived from the 
characteristics of the solar noise. The evaluation of B^ 1 ' 1 ' 
is less straightforward than that of C^ 1 ' 1 ) because we need 



3.3.4. Probability density of the observation and likelihood 

The statistical distribution of the Fourier spectra or of the 
vector z y is a multi-normal distribution. The probability 
density is given by: 



PvW) = 



e 2 



i4»v-V)*«M 



(31) 



where d is the number of elements of z y , V is a short no- 
tation for the following matrix: \v) = M^' ( \v) + 
B^ 1 ' 1 \v); this is the matrix given by the sum of thc p- 
mode and noise covariance matrix; the p modes and thc 
noises are assumed to be independent of each other. The 
matrix V^' \v) can also be built from sub-matrices as: 
y(M') = j^(i,0 + ,g(M') ; as a resu i t y(i,0 i s a l so hermi- 
tian. Equation (31) is the most general formulation for any 
multi-normal distribution with a given covariance matrix 
V (Kendall and Stuart, 1967). 
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Using Eq. (31), we can write the likelihood L of an 
observation of z y {vi) at N different frequencies vi as given 
by: 



I*™ 



n 



(32) 



We assumed that the frequency bins are independent of 
each other. This is the case when the data have no gaps. 
For unresolved observation having gap, the expression of 
the likelihood becomes extremely complicated as shown 
by Gabriel (1994) . For resolved observation having gaps, 
as for the LOWL data of Tomczyk et al. (1995), it is im- 
practicable to use the full formulation of the likelihood: 
Tomczyk et al. (1995) used Eq. (32) as an approximation 
for fitting the LOWL data. 

In principle, given the observed vector y, it is always 
possible in the absence of noise to recover the vector x. 
Due to the presence of noise only a solution close to the 
ideal one can be found that will minimize the correlation 
between the components. Provided that the leakage ma- 
trix can be inverted, we have by analogy to Eq. (12): 



(33) 



where C — C^ 1 ' 1 \ Then we can write a similar equation 
for z y and z$, as: 

z 5 = C^Zy (34) 



where C is defined as: 

r(M')f,/i - ( C ^ v ) C ^ v ) 
W \-Ci{y) CM 



(35) 



C^'' ) is a super matrix where C r and Ci are the real and 
imaginary parts of the complex matrix C"' 1 '. Using Eq. 
(34) to replace z y by z$ in Eq. (32) we can rewrite this 
latter as: 



N e -3*«("i)V'- 1 (fi)*«(i' i ) 1 



! £i M ' } (36) 



with V given by: 

V = C'W 1 = C- 1 mW)C t_1 + C^B^C 1 " -1 (37) 

We recognize in Eq (36) the probability density of the 
vector zg.{y) to a constant (i.e. |C| _iV ). As a matter of 
fact, it is well known that using a linear transformation 
similar to that of Eq. (34) will produce the new covariance 
matrix V of z% as written in Eq (37) (Davenport and 
Root, 1958). It can be easily shown using Eqs. (23) and 



(37) that the matrix D(v) = C^M^C 1 
and its element are given by: 

D m », m »(i/) = f m ,,(v) 



-l . 



is diagonal 



(38) 



where I" = I or V and m" = -I", . . . ,1" . Therefore Eq. 
(37) is the sum of a diagonal matrix representing the corre- 
lation between the p modes; and of a new noise covariance 
matrix representing the correlation of the components of 
the vector x after the transformation of Eq. (33). It means 
that x has no correlation due to the p modes as we could 
expect from Eq. (33): the leakage matrix of x is the iden- 
tity matrix. In summary, there is no gain in fitting data 
for which the leakage matrix is the identity matrix: the 2 
approaches are identical. The main problem is really to 
know the leakage matrices, not only theoretically but also 
experimentally: this is the subject of the Part II. 

It can be derived from Eq. (37) that it is also possible 
to remove correlation due to the noise by replacing C by 
a proper matrix associated with B^>' \ The derivation of 
this matrix is given in Appendix A. 



3.3.5. The use of the likelihood in practice 

When a single degree is observed, it is quite simple to 
maximize the likelihood of Eq. 32 using y, or using x as 
in Eq. 36. For low degree and low frequency modes, this 
is possible for I = 0,2,3. As soon as the mode linewidth 
increases, at high frequencies, the assumption of a single 
degree is not valid anymore. For example, I = and I = 1 
overlap with I = 2 and 1 = 3, respectively. At high fre- 
quencies, the effect of the aliasing degree should be taken 
into account. 

For the other low degree modes, the likelihood becomes 
somewhat more complicated. It is well known, that in the 
(to, v) diagramme of I = 1, there are leaks coming from 

1 = 6 and I = 9; in the (m, v) diagramme of I = 4, there 
are leaks of I = 7 and vice versa (Appourchaux et al. , 
1997). The leaks have severe effects on determination of 
the p-mode parameters of the 1 = 1. 

When many degrees are overlapping, one should use 
Eq. (32) using the covariance matrix for I and I'. Never- 
theless, we do not advice to do so for fitting the p modes; 
it has some severe computer speed penalty. Instead we ad- 
vice to clean the data by inverting the full leakage matrix 
taking into account the effects of the various degrees on 
each other, in a similar way to Eq. (33). This technique 
has been applied to the LOI and GONG data, and is de- 
veloped in Part II. 

Last but not least, when the signal-to-noise ratio is 
high (i.e. we neglect B^'' ) in Eq. 37), the elements of the 
vector Zx{y) are all independent of each other, leading 
to a statistical distribution which is a product of x 2 with 

2 degree of freedom. This is an approximation which is 
useful and less incorrect that using this statistics for the 
GONG data for the vector z y {v) as in Hill et al. 1996. 
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4. Monte-Carlo simulations 

4-1. Why are they needed? 

Before applying Eq. (32) to real data, it is always advisable 
to test the power of MLE on synthetic data, i.e. performing 
Monte-Carlo simulations. They are not merely for playing 
games; these simulations are real tools for understanding 
what we fit and how we fit it. Assuming that the statistics 
of the real solar spectra is known, performing Monte-Carlo 
is useful for the following reasons: 

— Assessing the model of the mode and noise covariance 

— Assessing the statistical distribution of the parameters 

— Assessing the precision of parameters 

First, the model of the covariances can be imperfect. The 
effect of an imperfect knowledge of the covariance can help 
us understand how these will influence the determination 
of the parameters, i.e. deriving the sensitivity of the sys- 
tematic errors to this imperfect knowledge. Second, the 
parameters derived by the MLE should have the desir- 
able properties of having a normal distribution; if not we 
advise to apply a change of fitted parameters. For exam- 
ple, as we will see later on, we do not fit the linewidth 
itself but the log of the linewidth. A normal distribution 
is necessary to derive meaningful error bars, this is the 
assumption behind Eq. (5). Third, in order to be able to 
derive a good estimate of the error bars using one real- 
ization, the standard deviation of a large sample of fitted 
parameters should be equal to the mean of formal errors 
return by the fit (See Eq. (4)). 

4-2. Generation of synthetic data for the LOI 

The performance of this instrument has been described in 
Appourchaux et al (1997). Briefly, it is a small instrument 
made of 12 pixels for detecting solar intensity fluctuations. 
The p-mode signals were generated in the Fourier spectra 
by using the following: 

N pix 

y{v)=CWx{v)+<Ty l i p i (39) 

i=l 

where y is the observed vector of 21 + 1 Fourier spec- 
tra, C^'^ is the leakage matrix given by Eq. (19), x is 
a complex random vector with 21 + 1 components (each 
component represents the signal of an I, m mode, with un- 
corrected real and imaginary part), y\ are computed as 
in Eq. (20) using spherical harmonics, and pi is the noise 
for a given pixel i. The variance of the real or imaginary 
part of the m-th component of x is given by f l m {v); the 
mean of x is 0. The function f l m {u) describes the profile 
of each m which is displaced from m by an amount which 
is given by: 

5 

Vm = vo + l^aiPf\m/l) (40) 



where the V\ are derived from the Clebsch-Gordan co- 
efficients, the expression of which can be found in Ritz- 
woller and Lavelly (1991); they are normalized such that 
Vf\l) = 1. Here we assumed a common linewidth for the 
l,n mode, and different amplitudes for the 21 + 1 com- 
ponents. The profile are symmetrical in the shape of a 
lorentzian. 

The variance of the pixel noise is assumed to be the 
same for the pixels with the same shape. The mean of the 
pixel noise is 0. For the LOI with its 12 pixels, there are 
3 different shapes giving 3 independent noises. 

After generating the synthetic signals according to Eq. 
(39), the data are fitted by minimizing the likelihood of 
Eq. (32). Figure 1 shows an example of Fourier spectra 
generated synthetically. The typical signal-to-noise ratio 
in the power spectra is about 20-30. The frequency res- 
olution is equivalent to 4 months of data. We performed 
1000 simulations of the spectra. 

4.3. Results 

4.3.1. For the nominal leakage matrix 

The data are fitted assuming a perfect knowledge of the 
leakage and noise covariance matrices, i.e. we know what 
we fit. Figure 2 shows the distribution of the fitted pa- 
rameters: the central frequency, splitting, log(linewidth), 
3x log(amplitude), 3x log(pixel noise). For the last 7 pa- 
rameters, we fit the log of the parameter because this 
transformation give a statistical distribution closer to a 
normal distribution (or log-normal distribution). It can 
be observed that the parameters derived are in most cases 
unbiased. Figure 3 shows the distribution of the error bars 
returned by the fit. In most case the mean of the error bars 
(returned by the fit) is not very different from the l-cr de- 
viation of the parameter distribution. Similar simulations 
have been performed for various degree (up 1 — 3). They 
show the same typical results as for Figs. 2 and 3, i.e. the 
fitted parameters are not, or weakly, biased, and the error 
bars returned by the fit give a good estimate of the real 
error bars. 

4.3.2. Influence of a wrong leakage matrix 

As was shown by Eq. (36), fitting p-mode spectra for 
which the leakage matrix is explicitly diagonal is equiv- 
alent to fitting p-mode spectra for which the matrix is not 
diagonal. Of course, it is always possible to construct data 
with a purely diagonal leakage matrix using Eq. (33), but 
we do so assuming that we know the leakage matrix C. As 
a matter of fact, what matters is not to have the identity 
matrix as leakage matrix, but more the knowledge of the 
latter. 

Hereafter, we have investigated the influence of a 
wrongly assumed leakage matrix on the fitted parameters 
of 1=1. We made 100 realizations and change the leakage 
parameter between m = — 1 and m = +1 by ±50% from a 
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M... A.. 



Frequency (in LiHz) 




Frequency (in //Hz 



Fig. 1. (Left) Power spectra of a synthetic Z = 1 as it would be observed by the LOI. The frequency resolution corresponds to 
4 months of data. The signal-to-noise ratio is about 20-30. The traces from bottom to top corresponds to m=-l,0,+l. (Right) 
Fourier spectra for Z = 1 (same data). The first, third and fifth traces from the bottom represents the real part of the spectrum 
of m = —1,0 and 1, respectively; the other traces are the imaginary parts. The leakage between m = — 1 and m = +1 is 0.45 
in the Fourier spectra. 
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Fig. 2. Histograms for the fitted parameters: (Plain line) Data, (Dashed line) Normal distribution with the same mean and 
a as the fitted parameters. (Top) Frequency (in ^iHz), splitting a\ (in /iHz), log(7) (7 in /iHz); (Middle) log( Amplitude) for 
m=-l,0,l; (Bottom) log(pixel noise). For each histogram, the target value, the mean fitted value and the 1-a fitted valued are 
displayed. The Kolmogorov-Smirnov test (Kol.) is displayed for each histogram; a number close to show that the distribution 
is not normal. 




Fig. 3. Histograms for the error bars: (Plain line) Data, (Dashed line) Normal distribution with the same mean and a. (Top) 
Frequency error (in /iHz), splitting ai error (in /iHz), log(7) error (7 in /iHz); (Middle) log (Amplitude) error for m=-l,0,l; 
(Bottom) log(pixel noise) error. For each histogram, the target value, the mean fitted value and the 1-a fitted valued are 
displayed. The Kolmogorov-Smirnov test is displayed for each histogram; a number close to show that the distribution is not 
normal. 



nominal value for the LOI of 0.45. Figure 4 shows the influ- 
ence of varying the assumed leakage element on the fitted 
parameters. It is quite interesting to note that the inferred 
central frequency is insensitive to mistakes in the leakage 
matrix. The linewidth becomes underestimated when the 
error is larger than 20 %, while the amplitudes become 
overestimated. The most important result is the fact that 
the systematic error made on the splitting is not linear 
but quadratic. This systematic error can become as large 
as the error bars. For example, with 1 year of LOI data 
and averaging over 10 modes, the error bars on the mean 
splitting is about 15 nHz; this should be compared to a 
systematic error of 10 nHz for an error of 10% of the I = 1 
leakage elements. 

Another test similar to that of the I = 1 was performed 
with the / = 2 mode. We have assumed that all the off- 
diagonal elements of the leakage matrix were wrong by 
the same fixed amount. Figure 5 shows the results only 
for the splitting coefficients (from a\ to (14). The other 



parameters linewidth, amplitudes and noises behave in the 
same manner as for I = 1. The systematic error on the 
splitting has also the same quadratic dependence as for 
1 = 1. For I = 2 the splitting error bars are typically \/5 
smaller than for / = 1. In this case the systematic errors 
become larger than the error bars, and therefore start to 
influence the inverted solar rotation. 

It means that it is quite easy to underestimate the 
splitting whenever we under- or overestimate the leakage 
element. As a matter of fact, this behaviour was also found 
in the GONG data for I = 1 and 2 (Rabello-Soares and 
Appourchaux, 1998, in preparation). On the other hand, 
errors in the leakage matrix will not result in overestimat- 
ing the splitting. If the splitting is overestimated, the most 
likely source should be the presence of other degrees not 
taken into account in the analysis. 

We also checked the correlation of the splitting coeffi- 
cients derived for 1=2 . Figure 6 and 7 show respectively 
the variance and the covariance of the splitting coefficients 
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as a function of the leakage elements error. It can be con- 
cluded that the splitting coefficients become correlated 
only when a large overestimation of about 50% is made 
for the off-diagonal leakage elements. This result is only 
valid when fitting Fourier spectra. For other methods, such 
as fitting power spectra, possible correlation amongst the 
splitting coefficients could have drastic consequences for 
the inverted solar rotation profiles. 

5. Conclusion 

We have given a step by step recipe for fitting (m, v) 
Fourier spectra. If one wants to implement similar fit- 
ting technique, one should compute, first the leakage 
matrices according to Eqs. (13) and (14), second the 
mode covariance matrices with Eq. (22) (or using Eq. 
(23)), third the noise covariance matrices with Eq. (27) 
using a model of solar noise, fourth compute the like- 
lihood function using Eq. (32). The use of Monte- 
Carlo simulations will ensure the success of the imple- 
mentation. Routines for fitting p-mode Fourier spectra 
are available as freeware on the VIRGO home page: 
ftp://ftp.estec.esa.nl/pub/loitenerife/html/software.html; 
they are written in the IDL language. 

Last but not least, Eq. (36) showed us the equivalence 
between fitting data for which the leakage matrix is not 
the identity, and fitting data for which it explicitly is. This 
last statement is true provided that we know perfectly well 
the leakage matrix. In this case the p-mode parameters fit- 
ted using MLE are not, or very weakly, biased, and have 
minimum variance. We have also studied the effect of an 
imperfect knowledge of the leakage matrix on the fitted 
parameters, in order to derive the effect of systematic er- 
rors on the most interesting parameters: the splitting and 
mode frequency We found that the central frequency is in- 
sensitive to systematic errors in the leakage matrix, while 
the splitting coefficients (dj) have a quadratic dependence 
upon those errors. These systematic errors will have influ- 
ence on the inverted solar rotation profiles. 

Finally, we would like to stress again that the correct 
statistical treatment of the p-mode data is of vital impor- 
tance for deducing unbiased p-mode parameters. 

Acknowledgements. Many thanks to Takashi Sekii for construc- 
tive comments on the manuscript, and for extensive cyberspace 
chatting :-) . 

A. Appendix 

The purpose of this appendix is to show that using a 
proper matrix Cb, the noise covariance matrix of Eq. (37) 
given by: 

B '(M') =Cb 1 B (M') C T-i (A1) 

can have a diagonal form. The matrix B^'*') can be diag- 
onalized and we can write. 

B (j,n = p-i b (M') P T- 1 (A2) 



where b^ ,z ' is diagonal and P is an orthogonal matrix 
(p-i = piy R e pi ac i n g Eq. (A2) into Eq. (Al), we have: 

B'('.0 = C - 1 p- 1 bW)p T " 1 Cr 1 (A3) 

Since B^'' ) is positive definite all its eigenvalues are pos- 
itive, therefore the square root of b^' Z ) is defined. There- 
fore if we apply the following transformation to the data: 

C B = p-VbW) (A4) 
we can rewrite Eq. (A3) as: 

B ,(M,) = I (A5) 

where I is the identity matrix. So replacing C in Eq. 33 by 
Cb will have the effect of removing the artificial correlation 
due to the noise, and also of performing a normalization. 
We should point out that the transformation matrix Cb 
that can achieve this is not unique, and any multiplication 
by an orthogonal matrix will achieve this. Nevertheless, we 
give a solution to the problem which can be solved as an 
eigenvalue and eigenvector problem. The transformation 
given above does not remove the artificial correlation due 
the p modes but more or less preserve it. This can have 
some useful application when one wants to produce spec- 
tra with uncorrelated noise but with correlated p-mode 
signals. 
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Fig. 4. Influence of the fitted parameters to relative changes of the assumed leakage element between m — — 1 and m = +1 
for I = 1. (Top) Frequency, splitting oi, log(7); (Middle) log (Amplitude) for m=-l,0,l; (Bottom) log(pixel noise). The target 
parameters are the same as for Fig. 2. Please note the parabolic shape for the splitting. 
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Fig. 5. Influence of the fitted splitting parameters to relative changes of the assumed of the assumed off-diagonal leakage 
element for I = 2. (Top, left) ai, target value: 410 nHz; (Top, right) ct2, target value: -30 nHz; (Bottom, left) (13, target value: 
-10 nHz; (Bottom, right) 04, target value: +50 nHz. 
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Fig. 6. Diagonal elements of the covariance matrix of the splitting coefficient, for I — 2. They are given as a function of the 
relative change of the assumed off-diagonal leakage element. (Top, left) For 01; (Top, right) For a 2 ; (Bottom, left) For a 3 ; 
(Bottom, right) For 04. 
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Fig. 7. Off-diagonal elements of the covariance matrix of the splitting coefficient, for / = 2. They are given as a function of the 
relative change of the assumed off-diagonal leakage element. (Top, left) For ai and a 2 ; (Top, right) For ai and a 3 ; (Middle, left) 
For ai and 04; (Middle, right) For 02 and ay, (Bottom, left) For 02 and 04; (Bottom, right) For 03 and 04. 



